function [sol] = Solveur_primal_GCP(nombreElements,nombreDomaines, F, donnees, errMax)

keq = donnees.E*donnees.S*nombreElements*nombreDomaines/donnees.L;

% creation des matrices blocs
Kii = keq*2*eye(nombreElements-1);
Kiidess = -keq*eye(nombreElements-2);
Kii(2:end,1:end-1) = Kii(2:end,1:end-1) + Kiidess;
Kii(1:end-1,2:end) = Kii(1:end-1,2:end) + Kiidess;

Kbb = keq*eye(2);

Kbi = zeros(2,nombreElements-1);
Kbi(end,end)= -keq;
Kbi(1,1)= -keq;

Kelem = donnees.E*donnees.S*nombreElements/donnees.L*[1,-1;-1,1];

Vtmp = Kii\Kbi';
Ktmp = Kbi*Vtmp;

% creation du complement de Schur primal elementaire
Spelem = Kbb-Ktmp; 

%% Construction du projecteur

% Calcul de G
    % Construction de A
    A = zeros(nombreDomaines,2*nombreDomaines-1) ;
    for i=1:nombreDomaines-1
        A(i,2*i-1) = 1 ;
        A(i,2*i) = 1 ;
    end
    A(end,end) = 1 ;
%     A = A(2:end,2:end) ;
    A 
    
    % Construction de A~
    Atilde = 0.5*A ; %Vrai dans le cas particulier de notre poutre en traction
    Atilde(end,end) = Atilde(end,end)*2 ;
    
    % Construction des modes de corps rigide
    R = zeros(2*nombreDomaines-1,nombreDomaines) ;
    for i=2:nombreDomaines
        R(2*i-2,i) = 1 ;
        R(2*i-1,i) = 1 ;
    end
    R(1,1) = 1 ;
    R
%     R = R(2:end,:) ;
    
    G = Atilde*R ;
    
    % Calcul du compélement de Schur assemblé
    Sp = Assemblage(Spelem,nombreDomaines) ;
    Sps = elimin(Sp) ;
    
    Gs = G
%     Gs = elimin(G) ;
    
% Q = G*inv(G'*Sp*G)*G' ;
% P = - Q*Sp ;
% P = P + eye(size(P(1))) ;

Q = Gs*inv(Gs'*Sps*Gs)*Gs'
P = - Q*Sps ;
P = P + eye(size(P,1)) ;

%  P = eye(size(P,1)) ;

%% Calcul de l'approximation de l'inverse de Sp

% Définition de t
t = zeros(2,nombreElements+1) ;
t(1,1) = 1 ;
t(end,end) = 1 ;
t

% Calcul du pseudo inverse de K
K = Assemblage(Kelem,nombreElements) ;
pinvK = pinv(K) ;

Ks = K(2:end,2:end) ;
pinvKs = pinv(Ks) ;

% Calcul du complement de Schur dual elementaire
Sdelem = t*pinvK*t' ;

% Assemblage du complement de Schur dual
% Sd = Assemblage(Sdelem,nombreDomaines) ;
% Sds = Sd(2:end,2:end) ;
Sd = Assemblage(Sdelem,nombreDomaines) ;
Sds = elimin(Sd) ;

% Calcul de l'approximationd de l'inverse de Sp
% iSp = 0.25*Sd ;
% iSps = iSp(2:end,2:end) ;
iSps = 0.25*Sds ;
iSps(end,end) = iSps(end,end)*4 ;

% initialisation

bp=[zeros(nombreDomaines-1,1);F];
u = Q*bp ;


niter = size(u,1) ;


r = P'*bp ;
z = iSps*r ;

w = zeros(size(z,1),niter+1) ;
p = zeros(size(z,1),niter+1) ;
alpha = zeros(niter+1,1) ;

w(:,1) = z ;

for i = 1:(niter)
    if norm(r,2)<errMax
        i
        break
    end
    
    p(:,i) = P'*Sps*w(:,i) ;
    alpha(i) = z'*r/(p(:,i)'*w(:,i)) ;
    u = u + alpha(i)*w(:,i) ;
    r = r - alpha(i)*p(:,i) ;
    z = iSps* r ;
    
    proj = zeros(size(z,1),1) ;
    for j = 1:i
        beta = - z'*p(:,i)/(p(:,i)'*w(:,i)) ;
        proj = proj + beta*w(:,j) ;
    end
    
    w(:,i+1) = z + proj ;
    
end

% reconstruction de la solution
u = [0;u] ;

sol= zeros(nombreElements*nombreDomaines+1,1);

for i = 1:nombreDomaines
    sol((i-1)*nombreElements+1)=u(i);
    Uinter = (-Kii\Kbi')*u(i:i+1);
    
    sol((i-1)*nombreElements+2:i*nombreElements)=Uinter;
    
end
sol(end) = u(end);

% r = zeros(nombreDomaines,1);
% 
% r(1) = bp(1)+Spelem(2,2)*u0(1);
% 
% for k=1:nombreDomaines-1
%     
%     r(k:k+1) =r(k:k+1) + bp(k:k+1) - Spelem*u0(k:k+1);
%     
% end
% 
% 
% w = r;
% 
% z = zeros(nombreDomaines,1);
% z(1) = Spelem(2,2)*w(1);
% 
% 
% 
% for k=1:nombreDomaines-1
%     
%     z(k:k+1) = z(k:k+1)+ Spelem*w(k:k+1);
%     
% end
% 
% 
% u = u0;
% B = 0;
% 
% % resolution par gradient conjugue
% 
% for i=1:nombreDomaines+1
%   a = dot(r,r)/dot(z,w); 
%   u = u+a*w;
%   d = dot(r,r);
%   r = r - a*z;
%     if( norm(r) < errMax )
%         break;
%     end
%     
%     B = dot(r,r)/d;
%     w = r + B*w;
% 
%     
%     z = zeros(nombreDomaines,1);
%     z(1) = Spelem(2,2)*w(1);
%     
%     for k=1:nombreDomaines-1
%         
%         z(k:k+1) = z(k:k+1)+Spelem*w(k:k+1);
%         
%     end
% 
% end
% 
% u = [0;u]
% 
% % reconstruction de la solution
% sol= zeros(nombreElements*nombreDomaines+1,1);
% 
% for i = 1:nombreDomaines
%     sol((i-1)*nombreElements+1)=u(i);
%     Uinter = (-Kii\Kbi')*u(i:i+1);
%     
%     sol((i-1)*nombreElements+2:i*nombreElements)=Uinter;
%     
% end
% sol(end) = u(end);
% 
% plot(sol);

end